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Abstract 


A  four-parameter  fit  is  developed  for  the  class  of  integrals  known  as  the  exponential  integral 
(real  branch).  Unlike  other  fits  that  are  piecewise  in  nature,  the  current  fit  to  the  exponential 
integral  is  valid  over  the  complete  domain  of  the  function  (compact)  and  is  everywhere  accurate 
to  within  ±0.0052%  when  evaluating  the  first  exponential  integral,  E^.  To  achieve  this  result, 
a  methodology  that  makes  use  of  analytically  known  limiting  behaviors  at  either  extreme  of  the 
domain  is  employed.  Because  the  fit  accurately  captures  limiting  behaviors  of  the  function, 
more  accuracy  is  retained  when  the  fit  is  used  as  part  of  the  scheme  to  evaluate  higher-order 
exponential  integrals,  E„,  as  compared  with  the  use  of  brute-force  fits  to  E^,  which  fail  to 
accurately  model  limiting  behaviors.  Furthermore,  because  the  fit  is  compact,  no  special 
accommodations  are  required  (as  in  the  case  of  spliced  piecewise  fits)  to  smooth  the  value,  slope, 
and  higher  derivatives  in  the  transition  region  between  two  piecewise  domains.  The  general 
methodology  employed  to  develop  this  fit  is  outlined,  since  it  may  be  used  for  other  problems 
as  well. 
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1.  Introduction 


A  number  of  useful  integrals  exist  for  which  no  exact  solutions  have  been  found.  In  other 
cases,  an  exact  solution,  if  found,  may  be  impractical  to  utilize  over  the  complete  domain  of  the 
function  because  of  precision  limitations  associated  with  what  usually  ends  up  as  a  series  solution 
to  the  challenging  integral.  For  many  of  these  integrals,  tabulated  values  may  be  pubhshed  in 
various  mathematical  handbooks  and  articles.  In  some  handbooks,  fits  (usually  piecewise)  also 
are  offered.  In  some  cases,  an  application  may  be  forced  to  resort  to  numerical  integration  in 
order  to  acquire  the  integrated  function.  In  this  context,  compact  {i.e.,  not  piecewise)  analytical 
fits  to  some  of  these  problematic  integrals,  accurate  to  within  a  small  fraction  of  the  numerically 
integrated  value,  serve  as  a  useful  tool  to  applications  requiring  the  results  of  the  integration, 
especially  when  the  integration  is  required  numerous  times  throughout  the  course  of  the 
application.  Furthermore,  the  ability  and  methodology  to  develop  intelligent  fits,  in  contrast  to 
the  more  traditional  “brate  force”  fits,  provide  the  means  to  minimize  parameters  and  maximize 
accuracy  when  tackling  some  of  these  difficult  functions.  The  exponential  integral  will  be  used 
as  an  opportunity  to  both  demonstrate  a  methodology  for  intelligent  fitting  as  well  as  for 
providing  an  accurate,  compact,  analytical  fit  to  the  exponential  integral. 

The  exponential  integral  is  a  useful  class  of  functions  that  arise  in  a  variety  of 
applications,  including  the  theory  of  interatomic  energy  potentials  [1],  the  theory  of  lethal  areas 
for  fragmenting  warheads  [2],  and  the  theory  of  ballistic  penetration  [3],  to  name  but  a  few.  The 
real  branch  of  the  family  of  exponential  integrals  may  be  defined  by 


Epi)  =  j£_  dt  ,  (1) 

X  ^ 


where  n,  a  positive  integer,  denotes  the  specific  member  of  the  exponential-integral  family.  The 
argument  of  the  exponential  integral,  rather  than  expressing  a  lower  limit  of  integration  as  in 
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eqn(l),  may  be  thought  of  as  describing  the  exponential  decay  constant,  as  given  in  this 
equivalent  (and  perhaps  more  popular)  definition  of  the  integral: 


(2) 


Integration  by  parts  permits  any  member  of  the  exponential-integral  family  to  be 
converted  to  an  adjacent  member  of  the  family,  by  way  of 


00 

^  ^  ) 

—  -  \—dt 

J  n 

X 

x”  J  t” 

K  J 

(3) 


expressable  in  terms  of  £„  as 


e"""  -  xEp:)  (n=  l,2, 3...) 


(4) 


Through  recursive  employment  of  eqn  (4),  all  the  members  of  the  exponential  integral  family  may 
be  analytically  related.  However,  this  technique  only  allows  for  the  transformation  of  one 
integral  into  another.  There  remains  the  problem  of  evaluating  E^(x).  There  is  an  exact  solution 
to  the  integral  of  appearing  in  a  number  of  mathematical  references  [4,5],  which  is 

obtainable  by  expanding  the  exponential  into  a  power  series  and  integrating  term  by  term.  That 
exact  solution,  which  is  convergent,  may  be  used  to  specify  E^{x)  as 


Ey{x)  =  -y  -  ln(x)  +  ^ 


2-2! 


3-3! 


(5) 
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Euler’s  constant,  y,  equal  to  0.57721...,  arises  in  eqn(5)  when  the  power-series  expansion  for 
{e~‘/t)  is  integrated  and  evaluated  at  its  upper  limit,  as  x  [6]. 

Employing  eqn  (5),  however,  to  evaluate  Ei(x)  is  problematic  for  finite  x  significantly 
larger  than  unity.  One  might  well  ask  of  the  need  to  evaluate  the  exponential  integral  for  large 
X,  since  the  function  to  be  integrated  drops  off  so  rapidly  that  the  integral  is  surely  a  very  flat 
function.  Such  reasoning  is  true  when  comparing  the  integrand  at  large  x  to  that  at  small  x. 
However,  the  definition  of  eqn(l)  has  as  its  upper  limit  not  a  small  value  of  x,  but  rather  that 
of  oo.  Therefore,  the  actual  values  for  E„(x)  are  extremely  small  numbers  for  large  values  of  x. 
Thus,  it  is  not  sufficient  merely  to  select  enough  terms  of  eqn  (5)  to  evaluate  the  integral  to 
within  a  value  of,  for  example,  ±0.0001  because  the  actual  integral  value  for  large  x  would  be 
smaller  than  this  arbitrary  tolerance.  To  draw  an  analogy,  it  would  be  like  saying  that  it  is  good 
enough  to  approximate  as  0.0  for  x>  10,  since  its  actual  value  is  within  0.0001  of  zero.  For 
some  applications,  such  an  approximation  may  be  warranted.  In  general,  though,  such  an 
approximation  is  mathematically  .unacceptable.  Worse  yet,  as  seen  from  eqns(l)  and  (2),  the 
need  to  evaluate  the  exponential  integral  for  large  arguments  can  arise  in  real-world  problems 
from  either  a  large  integration  limit  or  a  large  value  of  an  exponential  decay  constant.  Thus,  the 
need  to  evaluate  exponential  integrals  for  large  values  of  the  argument  is  established.  It  is  here 
that  the  practical  problems  with  the  evaluation  of  eqn  (5)  become  manifest. 

First,  the  number  of  terms,  N,  required  to  achieve  convergence  rises  rapidly  with 
increasing  x,  making  the  summation  an  inefficient  tool,  even  when  expressed  as  a  recursion 
relation  (for  three  digits  of  accuracy,  N  is  observed  to  vary  roughly  as  9  +  1.6x,  for  1  <x<7). 
More  important,  however,  is  the  fact  that,  for  calculations  of  finite  precision,  the  accuracy  of  the 
complete  summation  will  be  governed  by  the  individual  term  of  greatest  magnitude.  The  source 
of  the  problem  is  that  as  x  is  increased,  the  total  summation  decreases  in  magnitude  more  rapidly 
than  a  decaying  exponential,  while  at  the  same  time,  the  largest  individual  term  in  the  series  is 
observed  to  grow  rapidly  with  increasing  x  (~10'  for  x  =  7,  ~10^  for  x  =  10,  ~10^  for  x  =  13,  etc.). 
The  magnitude  of  this  largest  individual  term  consumes  the  available  precision  and,  as  a  result, 
leaves  little  or  none  left  for  the  ever-diminishing  net  sum  that  constitutes  the  desired  integral. 
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Literally,  the  use  of  eqn  (5),  even  with  (32-bit)  double  precision,  does  not  permit  the 
exponential  integral  to  be  evaluated  to  three  places  for  x:>  14  in  any  case,  and  with  the  situation 
worsening  for  lesser  precision.  For  these  reasons,  the  use  of  eqn  (5)  to  evaluate  the  exponential 
integral  numerically  for  large  x  is  wholly  unsuitable. 

Others  have  obviously  recognized  this  problem,  as  the  exponential  integral  in  some 
handbooks  [5]  receives  a  whole  chapter  of  attention.  For  large  x,  a  continued  fraction  exists  that 
converges  to  the  integral,  given  by 


E^{x)  =  e'^ 


1  + 


1  + 


X  + ... 


(6) 


But  as  X  becomes  smaller,  the  number  of  terms  required  for  convergence  rises  quickly.  Similar 
arguments  apply  for  the  use  of  an  asymptotic  expansion  for  E^,  which  also  converges  for  large 
X.  As  such,  the  more  typical  approach  employed  by  handbooks  is  that  of  a  fit.  While  some  steps 
are  taken  to  make  the  fits  intelligent  {e.g.,  transformation  of  variables),  the  fits  are  all  piecewise 
over  the  domain  of  the  integral. 

Cody  and  Thatcher  [7]  have  performed  what  is  perhaps  the  definitive  work,  with  the  use 
of  Chebyshev  approximations  to  the  exponential  integral  E^.  Like  others,  they  fit  the  integral 
over  a  piecewise  series  of  subdomains  (three  in  their  case)  and  provide  the  fitting  parameters 
necessary  to  evaluate  the  function  to  various  required  precisions,  down  to  relative  errors  of  10"^°. 
One  of  the  problems  with  piecewise  fitting  over  two  or  more  subdomains  is  that  functional  value 
and  derivatives  of  the  spliced  fits  will  not,  in  general,  match  at  the  domain  transition  point,  unless 
special  accomodations  are  made.  This  sort  of  discontinuity  in  functional  value  and/or  slope, 
curvature,  etc.,  may  cause  difficulties  for  some  numerical  algorithms  operating  upon  the  fitted 
function.  Numerical  splicing/smoothing  algorithms  aimed  at  eliminating  discontinuities  in  the 
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value  and/or  derivatives  of  a  piecewise  fit  are  not,  in  general,  computationally  insignificant. 
Problems  associated  with  piecewise  splicing  of  fits  may  also  be  obviated  by  obtaining  an  accurate 
enough  fit,  such  that  the  error  is  on  the  order  of  magnitude  of  the  limiting  machine  precision. 
This  alternative,  however,  requires  the  use  of  additional  fitting  parameters  to  acquire  the 
improved  precision.  Thus,  regardless  of  approach,  the  desire  to  eliminate  discontinuities  in  the 
function  and  its  derivatives,  between  piecewise  splices,  requires  extra  computational  effort.  One 
final  benefit  to  be  had  by  avoiding  the  use  of  piecewise  fits  is  the  concomitant  avoidance  of 
conditional  (Le.,  IF...THEN)  programming  statements  in  the  coding  of  the  routine.  The  use  of 
conditional  statements  can  preclude  maximum  computing  efficiency  on  certain  parallel  computing 
architectures. 

Therefore,  an  alternate  method  is  devised  to  approximate  the  exponential  integral.  A 
compact,  analytical  fit  that  strives  to  predict  the  exponential  integral  to  within  a  small  percentage 
of  the  actual  value  is  sought.  Ideally,  the  fit  will  remain  accurate  over  all  x,  which  will 
distinguish  it  fi'om  other,  piecewise,  fits  to  the  exponential  integral. 


2.  The  Approach 

Since  large  values  of  x  (x  >  1)  are  the  root  of  the  current  difficulty,  this  paper  first  focuses 
on  that  region  of  the  domain.  Fitting  any  function  accurately  over  the  expansive  domain  x>  1 
can  be  a  daunting  task.  Even  more  so  is  the  family  of  exponential  integrals,  eqn(l),  which 
experience  drastic  order  of  magnitude  changes  over  that  domain.  To  map  the  infinite  into  an 
acceptably  small,  finite  domain,  transform  the  exponential  integrals  by  way  of  the  substitution 
w=  1/x  (and  transform  the  dummy  variable,  t,  by  way  of  u=l/t).  In  this  manner,  the  domain 
x>  1  maps  into  the  transformed  domain  0  <  w <  1,  which  is  much  more  amenable  to  study.  (Of 
course,  the  downside  is  that  values  of  x<  1  transform  into  the  expansive  domain  w>  1.)  With 
this  simple  substitution,  the  family  of  exponential  integrals  given  by  eqn  (1)  becomes 

w 

E^{w  =  l/x)  =  — —  e'^'“  du  .  (7) 
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Eqn  (7)  is  then  numerically  integrated  with  respect  to  u  for  several  representative  values 
for  n  (using  Simpson’s  Rule,  with  a  step  size  of  1x10”’  across  the  small-w  domain),  so  as  to 
obtain  a  better  understanding  of  the  function.  At  small  w  (less  than  0.0016,  corresponding  to 
X  =  625),  the  ability  of  double  precision  (64  computer  bits)  to  represent  the  tiny  exponentials  is 
at  its  limit.  For  the  three  logarithmic  cycles  of  w  between  0.001  and  1,  the  exponential  integral 
varies  by  more  than  300  orders  of  magnitude.  Figure  1  shows  several  of  these  integrals  over 
40  of  those  orders  of  magnitude.  Clearly,  a  transformation  is  required  to  collapse  this  vast 
change  in  scale  into  something  less  dramatic.  Such  a  transformation  is  found,  when  it  is  realized 
that 

lim  =  1  .  (8) 

w  e 

Therefore,  E„  is  transformed  by  way  of  the  relation  F„  =  £„/(we"’'’").  In  this  way,  F„ 
equals  unity,  when  evaluated  at  w  =  0.  Figure  2  depicts  Fj,  F2,  and  F3  over  several  subdomains 
of  w.  It  is  observed  from  Figure  2a  that,  even  though  the  values  of  E„  cannot  be  evaluated  using 
double-precision  arithmetic  for  w  much  below  w  =  0.0016,  the  transformation  of  E„  into  F„  allows 
for  an  excellent  interpolation  between  the  (w,F„)  point  (0,1)  and  small  values  of  w  since,  unlike 
E„,  the  behavior  of  F„  at  low  w  is  nearly  linear.  As  observed  from  Figure  2a,  the  initial  slopes 
of  F„  Fj,  and  F,  appear  to  be  -1,  -2,  and  -3,  respectively.  To  obtain  these  values  exactly,  let 
F„=l-a„w  +b„ for  infinitesimal  w.  Substituting  the  definition  for  F„  (and  with  the  use 
of  eqn  [7])  results  in 

W 

E^  =  du  =  (1  -  aw  +  bw^  -  ...)  .  (9) 

0 


By  taking  the  derivative  of  this  expression  and  collecting  terms,  a„  must  equal  n,  and  b„  must 
equal  n{n+l)  in  the  small-w  limit.  These  derived  values  for  a„  confirm  the  slope  observations 
from  Figure  2a. 

Rather  than  attempting  a  fit  to  the  F„  directly,  however,  the  quantity  Ffw  is  first  plotted 
in  Figure  3.  At  large  w,  the  curve  appears  to  be  logarithmic  in  nature,  in  comparison  to  the  curve 
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Figure  1.  The  exponential  integrals  and  expressed  in  terms  of  h’  =  l/x,  in  relation  to 
the  decaying  exponential  function. 


Figure  2.  The  functions  F.,  F2,  and  7^3,  (a)  in  the  small- w  limit,  corresponding  to  x  — >  <»; 

(b)  over  the  domain  0  <  tv  <  1,  corresponding  to  x  >  1;  (c)  over  six  logarithmic 
cycles  of  >v,  corresponding  to  0.001  <x<  1000. 
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Figure  3.  A  comparison  of  Fj  •  w  to  ln(l+0.562»r),  for  larger  w,  above  unity. 
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ln(l  +  0.562  w),  shown  in  the  same  figure.  Large  values  of  w  correspond  to  small  x,  where 
eqn  (5)  is  more  easily  evaluated.  From  the  nature  of  eqn  (5)  as  x  approaches  zero,  it  is  seen  that 

lim  (x)  =  -y-  ln(Ac)  .  (10) 

x-^0 


From  eqn  (10),  by  transforming  x  to  1/w  and  to  a  comparison  of  the  term  F^'w 

to  the  functional  form  ln(l+^),  as  w— >oo,  reveals  that  k  must  take  on  the  exact  value 
0.5615... 


However,  for  small  w,  the  form  \n{\+e^w)  is  not  appropriate.  Rather,  recall  that,  for 
small  w,  the  functions  F„  are  shown  to  be  of  the  limiting  form  \-nw  +  n{n+\)w^.  Thus,  for  the 
case  of  n  =  1  (the  lowest-order  exponential  integral),  Fj  •  w  must  follow  as  w  -  for  small 

w.  A  Taylor-series  expansion  of  \n.{\-¥e^\v)  about  w  =  0,  by  contrast,  is  given  (to  two  terms)  by 
the  form  e^w-  1/2-  As  compared  to  w  -  +  2h^,  the  Taylor  expansion  is  observed  to 

have  the  wrong  initial  slope  (e^,  instead  of  1.0)  and  curvature  (— as  compared  with  the  proper 
value  of  -2)  in  the  low-w  limit.  To  remedy  the  slope  and  curvature  mismatch,  compare  in 
Figure  4  the  function  F,*w  to  a  different  logarithmic  form  (one  possessing  the  proper  initial 
slope  of  1  and  curvature  of  -2) — ^namely,  the  function  ln[l+ln(l+w)],  which  provides  an  excellent 
low-w  approximation  to  F^-w  by  matching  the  value,  slope,  and  curvature  of  F, •  w. 

,  Therefore,  the  limiting  behavior  of  Fj  •  w  is  known  for  both  w  — >  0  as  well  as  w  — >  oo.  As 
such,  the  fitting  strategy  for  Ffw  is  now  apparent:  a  form  which  follows  ln[l+ln(l+w)]  for  low 
w  and  transitions  to  ln(l+e'^w)  for  larger  w.  There  are  a  number  of  ways  in  which  this  transition 
might  be  accomphshed,  but,  after  some  study,  the  following  explicit  form  that  exhibits  the 
sought-after  characteristics  has  been  chosen: 

Fj  -w  »  ln{  1  +  w  -  [w  -  hi(l  +  w)  ]  '/(w)  }  ,  (11) 
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Figure  4. 


w 


A  comparison  of  F^-w  to  ln[l+ln(l+>r)]  and  ln(l+0.562M'),  for  small  ir,  below 
unity. 
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where  f(w)  is  the  function  to  parameterize  and  fit  One  alternative  strategy,  that  of  having  two 
fitting  functions  (/,  to  transition  down  the  low-w  form  and  to  activate  the  large-w  solution), 
is  discarded  for  two  reasons:  (1)  more  fitting  parameters  are  required  for  two  transition 
functions,  as  compared  to  one;  and  (2) /flnCl+w)  and/2-e“^w  interact  in  a  nonlinear  way, 
making  it  difficult  to  avoid  excessive  slope  and  curvature  oscillations  in  the  transition  region. 
By  contrast,  the  form  selected  in  eqn(ll)  has  three  primary  virtues  above  many  other  fitting 
forms:  (1)  it  requires  only  a  single  function  to  be  fit,  thereby  minimizing  the  number  of  fitting 
parameters;  (2)  as  w  increases,  the  logarithmic  term  in  brackets  naturally  exhibits  diminished 
importance  with  respect  to  the  linear  term,  thereby  providing  an  automatically  smooth  transition 
to  the  high-w  solution;  and  (3)  because  at  small  w  the  logarithmic  term  behaves  like  w  and 
monotonically  diverges  from  it  thereafter,  the  value  and  every  higher  derivative  of  the  quantity 
[w-ln(l+>v)]  exhibit  smooth,  monotonic  behavior.  Such  monotonicity  helps  to  ensure  that  the 
behavior  of /is  also  smooth  and  monotonic. 

Recalling  that  Ffw  is  identical  to  E^le~^''^,  the  first  three  exponential  integrals  may  be 
given,  in  terms  of/(w),  as 


7 

£  =  - i/M  «  -Infl +w- [w -ln(l+w)] -/(w)  }  ,  (12a) 

i  u 


E=l.  \e-^"'du  «  —  (w  -  In  {1  +  w  -  [w  -  ln(l  +w)]  -/(w) })  ,  (12b) 

W  J  w 


and 


E=—  \u  e'^^^du  =  1 -  w  +  In  { 1  +  w  -  [w  -  ln(l  +w)]  -/(w)  })  .  (12c) 

'  J  2w^ 
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Higher-order  exponential  integrals  may  be  determined,  after  a  similar  fashion,  with  the  use  of 
eqn  (4),  expressed  in  w. 

3.  TheFit(s) 

From  eqn  (1 1),/(0)  =  1,  in  order  to  have  the  proper  low-w  behavior.  At  larger  w,/ should 
approach  the  constant  value  of  (1  -  =  0.4385...  Furthermore,  recall  that  for  infinitesimal  w, 

F,-w  follows  the  form  w  -m^  +  2w^.  And  though  the  low-w  stencil,  ln[l+ln(l+w)],  matches  the 
limiting  value,  slope,  and  curvature  of  F,  •  w,  the  initial  slope  of  /  may  be  further  restricted  to 
force  a  match  to  the  third  derivative  of  F,  •  w,  as  well.  By  equating  this  limiting  cubic  form  to 
the  right  side  of  eqn  (11),  the  third  derivative  may  be  taken  (tedious  though  it  may  be)  to  show 
that  12  =  -3^'  +  7,  or 

fo'  =  -5/3  .  (13) 

The  numerically  integrated  values  for  Fj  are  used  to  generate  explicit  values  for  /,  as  shown  in 
Figure  5.  To  achieve  a  fit  to  this  data,  an  appropriate  functional  form  is  needed  to  approximate 
the  explicit  data  shown  in  Figure  5.  Two  forms  are  now  presented.  The  simplest,  with  two  free 
parameters,  is 

.  (14) 

(1  +(A  +  5/3)w+  Bw^  ) 


This  function  begins  at  unity,  for  w  =  0,  has  an  initial  slope  of  -5/3,  and  asymptotes  to  (1  —  e"^) 
for  large  w,  thus  satisfying  both  the  small-  and  large-w  asymptotic  requirements.  The  value  of 
B  relative  to  A  helps  to  determine  the  magnitude  of  the  transition  region. 

The  quality  of  this  fit  may  be  compared  to  the  numerically  integrated  values  and  be 
expressed  as  a  maximum  percent-deviation  from  the  integrated  value,  as  in 

%  Error  =  100  x  max[(FJ^,  -  EJ  /  .  (15) 
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Figure  5.  A  graph  of  the  /  function  to  be  fit,  showing  a  limiting  value  of  unity  as  h-  0, 
and  the  asymptote, /=  (1  -  e"^),  corresponding  to  large  w. 
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Over  the  complete  domain  of  the  function,  the  resulting  fit  (attempting  to  minimize  the 
percent  error  in  £j)  to  this  two-parameter  version  of/ is  given  as 


r/  \  1  +  5.874w  +  (1  -  e'^)'10.8w^  ^  *u-  \ 

f(w)  =  - - - -  ,  0<w’<oo  (£^1  withm  ±0.067%).  (16a) 

1  +(5.874 +  5/3)w+  10.8w2 


Restricting  the  fit  to  the  subdomain  0<  w<  1  (jc>  1),  as  do  many  of  the  existing  fits,  results  in 


r/  \  1  +  4.311w  +  (1  -  e‘^)*6.851w^  ^ 

f{w)  =  - ^ - - -  ,  0<w<l  (£i  withm  ±0.002%). (16b) 

1 +(4.311 +5/3)>v+  6.85  lw2 


A  four-parameter  form  is  also  developed  for/,  given  by 


,,  ,  1  +  4.054W  +  (1 -e-^)  *6.207^2 

f{w)  =  - 1 - 

1  +  (4.054 +  5/3)w  +  [1  +  0.5032w/(29.3+iv)2]  -6.2071^2 


,  0<w<  oo 


(El  within  ±0.005 11%).  (16c) 


Expressing  eqn(12a)  in  terms  of  x,  Ei  may  be  given  as 


eo 

E/x)  =  j—  dt  ~  *  In  { 1  +  1/x  -  [  1/x  -  ln(l  +  1/x)]  -/(x)  }  , 


(17) 


with  /(x)  given  as 


m  = 


x"  +  5.874x  +  (l -e-^)-10.8 


x2  + (5.874 +  5/3)x+  10.8 


,0<x<oo  (£j  within  ±0.067%),  (18a) 


15 


/w  = 


x^  +  43nx  +  il  -e-Tf) -6.851 
jc2  + (4.311 +5/3)a:+  6.851 


,  1  <jc<oo  (£■,  within  ±0.002%), 


(18b) 


or 


4.054;c  +  (1 -^-^) -6.207 

(4.054  +  5/3);c  +  [1  +  0.5032x/  (1  +29.3a:)''  ]  -6,207 


,  0<JC<  oo 


(El  within  ±0.005 11%).  (18c) 

Expressed  in  FORTRAN,  this  fit  to  £,  may  be  coded  in  a  single  statement  as 

El  =  EXP(-X)  *  LOG(l.  +  l./X  -  (l./X  -  LOG(l.  +  l./X)) 

&  *(X**2  +  4.054*X  +  2.72202  ) 

&  /(X**2  +5.72067*X  +(1.  +  0.5032*X/(1.  +  29 . 3*X) **2 ) *6 . 207 )  ) 

To  reiterate,  this  fit  is  vaUd  over  the  real  domain  of  j:>0,  to  within  an  error  of  ±0.00511%. 
There  is  no  need  for  conditional  branching  based  on  the  value  of  x.  Indeed,  there  is  hardly  a 
need  to  make  the  calculation  external  to  the  calling  routine. 

4.  Quality  of  Fit 

The  accuracy  for  each  of  the  fits  to /is  expressed  in  eqn(18)  in  terms  of  the  maximum 
percent  deviation  of  the  fitted  E,  from  its  numerically  integrated  value.  Because  of  the  manner 
in  which  E,  terms  interact  with  higher-order  exponential  integrals  [see  eqn(4)],  the  error 
produced  when  using  the  fitted /to  compute  these  higher-order  integrals  will  differ  from  that  of 
Ej.  Using  the  four-parameter  eqn(18c)  to  conduct  the  comparison  of  Ej,  E2,  and  E^  (i.e., 
eqns  [12])  with  the  integrated  data.  Figure  6a  reveals  errors  never  exceeding  ±0.15%.  Recall  that 
any  one  of  these  exponential  integrals  varies  by  more  than  300  orders  of  magnitude  over  the 
domain  shown  (0.001  <  w  <  1000)  in  Figure  6a.  Figure  6b  shows  the  same  data  (in  linear  scale) 
but  focuses  upon  the  domain  0<w<  1,  Over  this  domain,  corresponding  to  x>  1,  the  error  is 
strictly  less  than  ±0.005%  for  E^,  ±0.029%  for  E2,  and  ±0.15%  for  E^- 
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(b) 

Figure  6.  The  percent-error  of  the  compact  analytical  fit,  when  used  to  evaluate  E^,  E^,  and 
Ey,  (a)  shown  over  six  logarithmic  cycles  of  w;  (b)  over  the  domain  0  <  w  <  1. 
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By  contrast,  a  four-parameter  fit  to  F,  published  by  the  (former)  U.S.  National  Bureau  of 
Standards  (NBS)  [5]  (now  the  National  Institute  of  Standards  and  Technology),  given  by 

p,  _  +  2.334733a:  +  0.250621 

'  ~  +  3.330657a:  +  1.681534 

is  valid  only  for  the  limited  domain  of  w<  1  (x>  1)  and  has  an  absolute  error  band  in  Fj  of 
±0.00005  (expressed  in  terms  of  percent  error  in  F„  it  is  ±0.0051%).  The  four-parameter 
eqn  (18c)  betters  this  number  with  ±0.005%,  over  the  domain  a:  >  1,  as  shown  in  Figure  7.  Even 
the  two-parameter  eqn  (18b),  when  restricted  to  the  limited  domain  of  m'<  1  (a:>  1),  betters  the 
NBS  four-parameter  fit  by  a  factor  of  2.5.  To  achieve  a  better  accuracy  than  the  current  fit,  the 
National  Bureau  of  Standards  specifies  an  eight-parameter  fit  in  double  precision  which,  again, 
is  valid  over  only  a  limited  subdomain  of  the  function.  At  the  other  end  of  the  functional  domain 
(small  x),  compare  the  current  four-parameter  fit  to  eqn  (5)  when  truncated  to,  for  example,  the 
cubic  term.  In  this  case,  the  current  fit  is  more  accurate  for  a:  >0.2  (w<5). 

Furthermore,  because  the  NBS  fit  fails  to  achieve  the  proper  limiting  derivatives  as  w 
approaches  zero  (a:  <=°),  the  use  of  the  NBS  fit  to  acquire  the  higher-order  exponential  integrals, 

by  way  of  eqn  (4)  results  in  larger  errors  as  w  ->  0,  whereas  errors  with  the  current  fit  approach 
zero  as  w  0.  For  example,  E2  has  a  limiting  error  in  excess  of  -0.4%  with  the  NBS  fit 
(compared  with  a  maximum  error  of  ±0.029%  with  the  current  fit),  and  F3  has  a  limiting  error 
in  excess  of  +120%  (compared  with  a  maximum  error  of  ±0.15%  for  the  current  fit). 

A  comparison  of  the  current  fit  to  those  published  by  Cody  and  Thatcher  [7]  is  completely 
similar  to  that  of  the  NBS,  both  in  trend  and  in  order  of  magnitude.  To  exceed  the  accuracy  of 
the  current  fit,  six  parameters  are  required  for  the  subdomain  0  <  a:  <  1,  six  additional  parameters 
for  the  domain  1  <  a:  <  4,  and  4  more  parameters  for  the  domain  a:  >  4.  Furthermore,  each  of  those 
sixteen  parameters  is  specified  to  seven  or  more  digits  of  precision.  The  behavior  when  using 
the  fit  of  Cody  and  Thatcher  to  compute  higher-order  E„  is  also  similar  to  that  of  the  NBS; 
namely,  errors  become  much  larger,  more  quickly,  because  of  the  failure  of  Cody  and  Thatcher’s 
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Figure  7.  A  comparison  of  the  current  fit  to  the  four-parameter  NBS  fit,  over  the  limited 
domain  of  0  <  ir  <  1. 


19 


fit  to  as)rmptote  to  the  proper  limiting  derivatives.  This  particular  shortcoming  of  their  fit  can 
be  circumvented,  but  only  through  the  use  of  significantly  greater  number  of  parameters  in  the 
fit  to  the  jc  >  4  subdomain, 

5.  Fitting  Strategies 

A  few  words  may  be  in  order  here,  on  the  formulation  of  fitting  strategies  for  functions 
in  general.  With  enough  parameters,  any  function  may  be  fit  to  a  high  degree  of  accuracy.  The 
goal,  however,  is  to  achieve  this  fit,  while  minimizing  the  number  of  parameters  and 
computational  expense  while,  at  the  same  time,  maximizing  the  accuracy  and  the  domain  of 
applicability.  Depending  on  how  little  or  great  an  effort  is  warranted  for  the  problem  at  hand, 
some  or  all  of  the  following  suggestions  may  be  of  utility  in  formulating  an  intelligent  fit  to  a 
specified  function,  G,  numerical  values  for  which  are  known  only  in  tabular  form: 

1)  Compress  the  function’s  domain  and  range  of  interest  to  a  manageable  size  (i.e.,  analyze  and 
fit  a  transformed  function,  H,  if  easier  to  digest  than  the  original  function,  G). 

2)  Obtain  limiting  behaviors  of  the  original  (G)  or  transformed  (H)  function,  if  possible.  This 
includes  not  only  functional  values,  but  slopes  and  higher  derivatives.  Such  limits  are  needed 
if  one  wishes  the  fit  to  accurately  converge  in  the  limits  of  the  domain. 

3)  Convert  limiting  behaviors  into  functional  stencils  (e.g.,  hf^rge-^,  gasymptote^  Often, 

the  more  derivatives  that  the  stencil  can  match,  the  easier  it  will  be  to  subsequently  fit  the 
transition  function. 

4)  Develop  a  functional  form  for  G  or,  if  easier,  for  H  (in  terms  of,  hopefully,  a  single  transition 

function,  /)  that  is  able  to  effectively  transition  between  the  limiting  behaviors  of  the  function. 
At  this  point,  the  unknowns  are  not  individual  fitted  parameters,  but  instead  the  transition 
function(s).  A  simple  example  would  be  H  From  this  developed  form. 


20 


quite  a  bit  should  be  known  about  the  general  behavior  of  the  transition  function,/.  For  example, 
it  might  be  known  that/(0)  is  1,  with  an  initial  slope  of  -1,  and  which  asymptotes  at  large  x  to  0. 

5)  Explicitly  solve  for  tabulated  values  of  the  transition  function,  using  the  numerically  integrated 

or  tabulated  values  for  G  or  H,  and  the  calculated  values  for  the  stencils,  h^„aii-x>  Using  the 
above  example,/ =  The  goodness  of  a  particular  functional  form  for 

G  or  //,  in  terms  off,  may  be  assessed  by  observing  how  smoothly  the  calculated  values  of  the 
/function  behave.  For  example,  it  is  generally  easier  to  fit  a  monotonic  transition  function  than 
one  which  oscillates.  If  the  computed  values  for  /don’t  seem  to  follow  a  desired  (easily  fitted) 
form,  return  to  step  4  and  reformulate  a  new  functional  form. 

6)  Based  on  observations  of  the  behavior  of  the  transition  function,  /,  specify  a  specific  form  for 
the  transition  function  in  terms  of  parameters  to  be  fit.  For  example, /=A/x+Bx^.  The  observed 
behavior  of  the  tabulated  values  for  /  should  help  in  the  specification  of  a  good  parameterized 
form. 

7)  Fit  the  parameters  associated  with  the  transition  function.  Ample  numbers  of  coded  programs 
exist  to  assist  in  this  task,  or  else  a  simple  (if  inefficient)  fitting  routine  may  be  written  to 
minimize  a  specified  error  function  (e.g.,  least  squares). 

8)  Once  the  fitted  parameters  are  obtained,  evaluate  the  accuracy  of  the  fit.  If  quality  of  fitted 
results  are  below  expectations,  return  to  step  6  or  4  to  reformulate  new  trial  functions. 

In  the  current  work  of  fitting  step  1  is  accomplished  by  transforming  both  the  domain 
(from  X  to  w  in  eqn  [7])  and  the  range  (from  Ej  to  Fj  and  eventually  to  Fj-w).  In  step  2,  eqn  (9) 
is  employed  to  obtain  the  slope  and  derivatives  of  F  for  infinitesimal  w  and  eqn  (10)  is  used  to 
express  the  analytically  known  behavior  of  Fj  for  small  x  in  terms  of  large  w.  Step  3  is 
accomplished  by  formulating  the  stencils  ln[l+ln(l+w)]  for  small  w  and  ln(l+e'^w)  for  large  w. 
Eqn  (11)  constitutes  the  functional  form  developed,  as  part  of  step  4.  The  transition  function,/ 
is  explicitly  calculated,  per  step  5,  and  plotted  in  Figure  5  to  reveal  its  desirable  (easily  fitted) 
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characteristics.  The  transition  function  is  specified,  in  terms  of  parameters,  in  eqn(14),  in 
accordance  with  step  6.  Using  a  computational  routine,  the  parameters  to  the  transition  function 
are  fitted,  as  given  in  eqns  (16),  per  step  7.  The  accuracy  is  evaluated  and  found  acceptable,  thus 
completing  step  8  and  the  process. 

This  methodology  may  also  be  used  to  generate  fits  to  experimental  data,  if  something 
is  analytically  known  about  the  physical  process  which  created  the  data.  In  contrast  to  the 
methodology  described  here,  polynomial  fits  are  very  popular  (especially  in  the  experimental 
community)  because  of  the  trivial  cost  in  obtaining  the  fit.  Unfortunately,  polynomial  fits  are 
virtually  guaranteed  to  offer  zero  extrapolative  capacity  beyond  the  domain  over  which  they  were 
fit.  Furthermore,  they  offer  zero  insight  into  the  physical  behavior  of  the  unknown  function  or 
data,  and  generally  require  many  more  parameters  to  achieve  the  same  level  of  accuracy  as  an 
intelligent  fit,  as  described  here.  And  without  some  forethought,  functions  such  as  the 
exponential  integrals  are  very  costly  to  fit  directly  with  polynomials  because  of  their  exponential 
nature  [e.g.,  a  truncated  eqn  (5)  is  essentially  a  very  costly  polynomial  fit  to  the  exponential 
integral]. 

In  contrast  to  the  polynomial  fit,  the  process  described  here  offers  a  logical  methodology 
for  generating  fits  to  difficult  functions  or  data,  such  that  good  results  may  be  obtained  with 
fewer  parameters  and  at  a  reasonable  cost.  For  an  oft-used  function,  the  extra  effort  required  to 
compose  an  intelligent  fit  might  very  well  be  justified.  A  quote  attributed  to  Cauchy,  the  19th 
century  mathematician,  contends  that  “Give  me  five  free  parameters  and  I  will  give  you  the 
equations  for  an  elephant,  but  give  me  a  sixth  free  parameter  and  I  will  make  the  elephant  wag 
its  tail.”  Such  versatility  with  so  few  parameters  is  surely  possible  only  by  way  of  an  intelligent 
fit. 


6.  Summary 

A  compact  (not  piecewise),  analytical  fit  has  been  developed  for  the  class  of  integrals 
known  as  exponential  integrals.  The  fit  has  been  shown  to  match  the  numerically  integrated 


22 


exponential  integral,  £„  to  within  an  error  of  0.0052%,  over  the  complete  domain  of  the  function. 
Related  integrals,  E2  and  £3,  were  computed  to  within  0.029%  and  0.15%,  respectively,  when 
using  the  current  £j  fit  as  part  of  the  general  E„  solution.  Achieving  such  a  feat  is  complicated 
by  the  fact  that  the  value  of  these  integrals,  being  exponential  in  nature,  can  change  several- 
hundred  orders  of  magnimde  over  a  reasonably  small  section  of  the  integral’s  domain. 

By  avoiding  a  series  of  piecewise  fits  to  describe  the  complete  domain  of  the  function, 
it  is  guaranteed  that  the  current  fit  is  continuous  in  value,  slope,  and  all  higher  derivatives.  By 
contrast,  a  series  of  piecewise  fits  must  either  contend  with  discontinuous  values,  slopes,  and/or 
higher  derivatives,  or  else  employ  a  smoothing  algorithm  to  splice  the  piecewise  fits.  Such 
smoothing  algorithms  could  be  more  computationally  expensive  than  the  fit  itself  and,  therefore, 
need  to  be  considered  as  part  of  the  computational  burden.  Additionally,  the  avoidance  of 
piecewise  fitting  permits  the  coding  of  the  fit  into  a  programming  language  without  the  use  of 
conditional  {i.e.,  IF...THEN)  statements.  For  some  parallel-computing  architectures,  the  use  of 
conditional  statements  (required  of  piecewise  fits)  precludes  maximum  computing  efficiency. 

The  ability  to  achieve  a  good  fit,  encompassing  the  complete  domain  of  the  function,  with 
the  use  of  only  four  parameters,  stems  completely  from  the  fact  that  so  much  could  be 
analytically  deduced  about  the  behavior  of  the  function  in  large-x  and  small-x  limits,  prior  to 
acmally  composing  the  form  of  the  fit.  In  the  small-x  limit,  the  exact  limiting  form  of  the  result 
is  analytically  known.  In  the  large-x  limit,  a  functional  stencil  that  matches  value,  slope,  and 
curvature  of  the  acmal  exponential  integral  is  employed.  Furthermore,  this  stencil  for  the  large-x 
(small-w)  limit,  given  by  ln[l+ln(l+>v)],  does  not  numerically  explode  as  w  increases  but,  instead, 
smoothly  fades  as  a  diminishing  percentage  of  the  large-w  (small-x)  solution.  In  this  manner, 
special  accommodations  (such  as  piecewise  fitting  of  the  solution)  do  not  have  to  be  made  in 
order  to  transition  from  the  small-w  to  large-w  solution. 

The  accuracy  of  the  fit  exceeds  a  comparable  {i.e.,  four-parameter)  fit  published  by  the 
(former)  U.S.  National  Bureau  of  Standards  [5].  Furthermore,  the  NBS  fit  is  accurate  over  only 
a  specified  subset  of  the  functional  domain.  If  the  current  fit  is  restricted  to  the  same  functional 
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subdomain  as  the  NBS  fit,  only  two  parameters  are  required  to  significantly  exceed  the  accuracy 
of  the  NBS  fit.  And  because  the  NBS  fit  fails  to  properly  match  the  limiting  slope  and  higher- 
order  derivatives  of  the  exponential  integral,  £„  the  use  of  the  NBS  fit  to  evaluate  the  higher- 
order  E„  produces  a  greater  limiting  error  than  the  current  fit:  over  120%  in  the  case  of  £3,  as 
compared  with  0.15%  by  using  the  current  fit.  Similar  results  arise  when  comparing  the  current 
fit  to  the  four-parameter  fits  of  Cody  and  Thatcher  [7],  which  are  more  recent  than  the  NBS  fits. 
However,  if  the  number  of  parameters  and  piecewise  nature  of  the  fits  are  not  issues,  Cody  and 
Thatcher  still  provide  the  most  accurate  of  all  fits.  To  achieve  errors  on  the  order  of  one  part 
in  10^°,  though,  a  total  of  52  high-precision  fitting  parameters  are  required  on  their  part. 

The  methodology  used  to  generate  the  current  fit  was  summarized  in  a  loosely  codified 
form.  This  methodology  can  be  of  assistance  in  generating  fits  to  other  difficult  functions  or 
even  experimental  data,  if  something  is  analytically  known  about  the  physical  processes  that 
created  the  data. 
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2  AEROJET 

ATTN  J  CARLEONE 
S  KEY 

PO  BOX  13222 

SACRAMENTO  CA  95813-6000 

2  AEROJET  ORDNANCE 
ATTN  P  WOLF 
G  PADGETT 
1100  BULLOCH  BLVD 
SOCORRO  NM  87801 


2  ALLIANT  TECHSYSTEMS  INC 
ATTN  R  STRYK 
G  R  JOHNSON  MNl  1-1614 
600  SECOND  ST  NE 
HOPKINS  MN  55343 

1  MARVIN  L  ALME 

2180  LOMA  LINDA  DR 
LOS  ALAMOS  NM  87544-2769 

1  APPLIED  RESEARCH  ASSOC  INC 
ATTN  JEROME  D  YATTEAU 

5941  S  MIDDLEFIELD  RD  SUITE  100 
LITTLETON  CO  80123 

2  APPLIED  RESEARCH  ASSOC  INC 
ATTN  DENNIS  GRADY 
FRANK  MAESTAS 

SUITE  A220 

4300  SAN  MATEO  BLVD  NE 
ALBUQUERQUE  NM  87110 

1  ATA  ASSOCIATES 
ATTN  W  ISBELL 
PO  BOX  6570 

SANTA  BARBARA  CA  93111 

1  BATTELLE 

ATTN  ROBER  M  DUGAS 

7501  S  MEMORIAL  PKWY  SUITE  101 

HUNTSVILLE  AL  35802-2258 

3  BOEING  AEROSPACE  CO 

SHOCK  PHYSICS  &  APPLIED  MATH 
ENGINEERING  TECHNOLOGY 
ATTN  R  HELZER 
T  MURRAY 
J  SHRADER 
PO  BOX  3999 
SEATTLE  WA  98124 

1  BOEING  HOUSTON  SPACE  STN 
ATTN  RUSSELL  F  GRAVES 
BOX  58747 
HOUSTON  TX  77258 

1  BRIGS  CO 

ATTN  JOSEPH  E  BACKOFEN 
2668  PETERSBOROUGH  ST 
HERNDON  VA  20171-2443 
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1  CALIFORNIA  RSCH  &  TECHNOLOGY 
ATTN  M  MAJERUS 
PO  BOX  2229 
PRINCETON  NJ  08543 

1  CENTURY  DYNAMICS  INC 
ATTN  N  BIRNBAUM 
2333  SAN  RAMON  VALLEY  BLVD 
SAN  RAMON  CA  94583-1613 

1  COMPUTATIONAL  MECHANICS 
CONSULTANTS 
ATTN  J  A  ZUKAS 
PO  BOX  11314 
BALTIMORE  MD  21239-0314 

1  CYPRESS  INTERNATIONAL 
ATTN  A  CAPONECCHI 
1201  E  ABINGDON  DR 
ALEXANDRIA  VA  22314 

1  DEFENSE  TECHNOLOGY  INTL.  INC 
ATTN  D  E  AYER 
THE  STARK  HOUSE 
22  CONCORD  ST 
NASHUA  NH  03060 

1  DESKIN  RESEARCH  GROUP  INC 
ATTN  EDWARD  COLLINS 
2270  AGNEW  RD 
SANTA  CLARA  CA  95054 

3  DOW  CHEMICAL  INC 
ORDNANCE  SYSTEMS 
ATTN  C  HANEY 
A  HART 
B  RAFANIELLO 
800  BUILDING 
MIDLAND  MI  48667 

1  G  E  DUVALL 

5814  NE  82ND  COURT 
VANCOUVER  WA  98662-5944 

3  DYNA  EAST  CORP 

ATTN  P  C  CHOU 
R  CICCARELLI 
WFLIS 

3620  HORIZON  DRIVE 
KING  OF  PRUSSIA  PA  19406 


3  DYNASEN 

ATTN  JACQUES  CHAREST 
MICHAEL  CHAREST 
MARTIN  LILLY 
20  ARNOLD  PL 
GOLETA  CA  93117 

1  R  J  EICHELBERGER 
409  W  CATHERINE  ST 
BEL  AIR  MD  21014-3613 

1  ELORET  INSTITUTE 

ATTN  DAVID  W  BOGDANOFF  MS  230  2 
NASA  AMES  RESEARCH  CENTER 
MOFFETT  FIELD  CA  94035 

3  ENIG  ASSOaATES  INC 
ATTN  J  ENIG 
D  J  PASTIME 
M  COWPERTHWAITE 
SUITE  500 

11120  NEW  HAMPSHIRE  AVE 
SILVER  SPRING  MD  20904-2633 

1  EXPLOSIVE  TECHNOLOGY 
ATTN  M  L  KNAEBEL 
PO  BOX  KK 
FAIRFIELD  CA  94533 

1  GB  TECH  LOCKHEED 
ATTN  JAY  LAUGHMAN 
2200  SPACE  PARK  SUITE  400 
HOUSTON  TX  77258 

2  GB  TECH  LOCKHEED 

ATTN  LUCILLE  BORREGO  C23C 
JOE  FALCON  JR  C23C 
2400  NASA  ROAD  1 
HOUSTON  TX  77058 

6  GDLS 

38500  MOUND  RD 
ATTN  W  BURKE  MZ436-21-24 
G  CAMPBELL  MZ436-30-44 
D  DEBUSSCHER  MZ436-20-29 
J  ERIDON  MZ436-21-24 
W  HERMAN  MZ  435-01-24 
S  PENTESCU  MZ436-21-24 
STERLING  HTS  MI  48310-3200 
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2  GENERAL  RESEARCH  CORP 
ATTN  A  CHARTERS 
TMENNA 
PO  BOX  6770 

SANTA  BARBARA  CA  93160-6770 

1  GRC  INTERNATIONAL 

ATTN  TIMOTHY  M  CUNNINGHAM 
5383  HOLLISTER  AVE 
SANTA  BARBARA  CA  93 1 1 1 

7  INST  OF  ADVANCED  TECHNOLOGY 
UNIVERSITY  OF  TX  AUSTIN 
ATTN  S  J  BLESS 
J  CAZAMIAS 
J  DAVIS 
H  D  FAIR 
TMKIEHNE 
D  LITTLEFIELD 
M  NORMANDIA 
4030-2  W  BRAKER  LN 
AUSTIN  TX  78759 

1  INTERNATIONAL  RESEARCH  ASSOC 
ATTN  D  L  ORPHAL 
4450  BLACK  AVE 
PLEASANTON  CA  94566 

1  INTERPLAY 

ATTN  F  E  WALKER 
18  SHADOW  OAK  RD 
DANVILLE  CA  94526 

1  ITT  SCIENCES  AND  SYSTEMS 
ATTN  J  WELBECK 
600  BLVD.  SOUTH,  SUITE  208 
HUNTSVILLE  AL  35802 

1  R  JAMESON 
624  ROWE  DR 
ABERDEEN  MD  21001 

1  KAMAN  SCIENCES  CORP 
ATTN  DENNIS  L  JONES 
2560  HUNTINGTON  AVE  SUITE  200 
ALEXANDRIA  VA  22303 


7  KAMAN  SCIENCES  CORP 
ATTN  J  ELDER 
RICHARD  P  HENDERSON 
DAVID  A  PYLES 
FRANK  R  SAVAGE 
JAMES  A  SUMMERS 
TIMOTHY  W  MOORE 
THY  YEM 

600  BLVD  S  SUITE  208 
HUNTSVILLE  AL  35802 

3  KAMAN  SCIENCES  CORP 
ATTN  SHELDON  JONES 
GARY  L  PADEREWSKI 
ROBERT  G  PONZINI 

1500  GRDN  OF  THE  GODS  RD 
COLORADO  SPRINGS  CO  80907 

4  KAMAN  SCIENCES  CORP 
ATTN  NASrr  ARI 
STEVE  R  DIEHL 
WILLIAM  DOANE 
VERNON  M  SMITH 

PO  BOX  7463 

COLORADO  SPRINGS  CO  80933-7463 

1  D  R  KENNEDY  &  ASSOC  INC 
ATTN  D  KENNEDY 
PO  BOX  4003 

MOUNTAIN  VIEW  CA  94040 

1  KERLEY  PUBLISHING  SERVICES 
ATTN  G  I  KERLEY 

PO  BOX  13835 

ALBUQUERQUE  NM  87192-3835 

2  KTECH  CORPORATION 
ATTN  FRANK  W  DAVIES 
LARRY  M  LEE 

901  PENNSYLVANIA  NE 
ALBUQUERQUE  NM  87110 

1  LIVERMORE  SOFTWARE  TECH  CORP 
ATTN  J  O  HALLQUIST 
2876  WAVERLY  WAY 
LIVERMORE  CA  94550 
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I  LOCKHEED  MARTIN  MISSLE  &  SPACE 
ATTN  WILLIAM  R  EBERLE 
PO  BOX  070017 
HUNTSVILLE  AL  35807 

3  LOCKHEED  MARTIN  MISSILE  &  SPACE 
ATTN  M  A  LEVIN  ORG  81  06  BLDG  598 
M  R  MCHENRY 

T  A  NGO  ORG  81  10  BLDG  157 
111  LOCKHEED  WAY 
SUNNYVALE  CA  94088 

4  LOCKHEED  MISSILE  &  SPACE  CO 
ATTN  JOHN  R  ANDERSON 
WILLIAM  C  KNUDSON 

S  KUSUMI  0  81  11  BLDG  157 
J  PHILLIPS  0  54  50 
PO  BOX  3504 
SUNNYVALE  CA  94088 

1  LOCKHEED  MISSILE  &  SPACE  CO 
ATTN  R  HOFFMAN 
SANTA  CRUZ  FACILITY 
EMPIRE  GRADE  RD 
SANTA  CRUZ  CA  95060 

1  LOCKHEED  NASA  JSC 
SPACE  SCIENCE  BRANCH 
ATTN  JAMES  HYDE 
BOX  58561  MC  B22 
HOUSTON  TX  77258 

1  MCDONNELL  DOUGLAS 
ASTRONAUTICS  CO 
ATTN  B  L  COOPER 
5301  BOLSA  AVE 
HUNTINGTON  BEACH  CA  92647 

2  ORLANDO  TECHNOLOGY  INC 
ATTN  DANIEL  A  MATUSKA 
MICHAEL  GUNGER 

PO  BOX  855 

SHALIMAR  FL  32579-0855 

1  PHYSICAL  SCIENCES  INC 
ATTN  PETER  NEBOLSINE 
20  NEW  ENGLAND  BUS  CTR 
ANDOVER  MA  01810 


NO.  OF 

COPIES  ORGANIZATION 

3  PHYSICS  INTERNATIONAL 
ATTN  R  FUNSTON 
G  FRAZIER 
L  GARNETT 
PO  BOX  5010 
SAN  LEANDRO  CA  94577 

1  PRC  INC 

ATTN  J  ADAMS 

5166  POTOMAC  DR  #103 

KING  GEORGE  VA  22485-5824 

1  RAYTHEON  ELECTRONIC  SYSTEMS 
ATTNRKARPP 
50  APPLE  HILL  DRIVE 
TEWKSBURY  MA  01876 

1  ROCKWELL  INTERNATIONAL 
ROCKETDYNE  DIVISION 
ATTN  H  LEIFER 
16557  PARK  LN  CIRCLE 
LOS  ANGELES  CA  90049 

1  ROCKWELL  MISSILE  SYS  DIV 
ATTN  T  NEUHART 
1800  SATELLITE  BLVD 
DULUTH  GA  30136 

1  SAIC 

ATTN  MICHAEL  W  MCKAY 
10260  CAMPUS  POINT  DR 
SAN  DIEGO  CA  92121 

1  SHOCK  TRANSIENTS  INC 
ATTN  DAVE)  DAVISON 
BOX  5357 

HOPKINS  MN  55343 

2  SIMULATION  &  ENG  CO  INC 
ATTN  ELSA  I  MULLINS 
STEVEN  E  MULLINS 

8840  HWY  20  SUITE  200  N 
MADISON  AL  35758 

2  SOUTHERN  RESEARCH  INSTITUTE 
ATTN  LINDSEY  A  DECKARD 
DONALD  P  SEGERS 
PO  BOX  55305 

BIRMINGHAM  AL  35255-5305 
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5  SRI  INTERNATIONAL 
ATTN  JAMES  D  COLTON 
D CURRAN 
R  KLOOP 
R  L  SEAMAN 
D  A  SHOCKEY 
333  RAVENSWOOD  AVE 
MENLO  PARK  CA  94025 

2  TELEDYNE  BROWN  ENGR 
ATTN  EM  W  BOOTH 
MARTIN  B  RICHARDSON 
PO  BOX  070007  MS  50 
HUNTSVILLE  AL  35807-7007 

1  ZERNOW  TECHNICAL  SVCS  INC 
ATTN  LOUIS  ZERNOW 
425  W  BONITA  AVE  SUITE  208 
SAN  DIMAS  CA  91773 

1  SUNY  STONEYBROOK 

DEPT  APPL.  MATH  &  STAT. 
ATTN  J  GLIMM 
STONEYBROOK  NY  11794 


ABERDEEN  PROVING  GROUND 

9  DIR,  USARL 

ATTN:  AMSRL-WM,  I  MAY 
AMSRL-WM-BC,  A  ZIELINSKI 
AMSRL-WM-BD, 

R  PESCE-RODRIGUEZ 
A  KOTLAR 

AMSRL-WM-BE,  S  HOWARD 
AMSRL-WM-MB,  G  GAZONAS 
AMSRL-WM-MC,  J  M  WELLS 
AMSRL-WM-T, 

W  F  MORRISON 
TW  WRIGHT 


43  DIR,  USARL 

ATTN:  AMSRL-WM-TA, 

W  GILLICH 
S  BILYK 
M  BURKINS 
W  BRUCHEY 
J  DEHN 
GFILBEY 
W A  GOOCH 
HW  MEYER 
E  J  RAPACKI 
J  RUNYEON 
AMSRL-WM-TB, 

RFREY 
P BAKER 
R  LOTTERO 
J  STARKENBERG 
AMSRL-WM-TC, 

W  S  DE  ROSSET 
T  W  BJERKE 
R  COATES 
F GRACE 
K  KIMSEY 
M  LAMPSON 
D  SCHEFFLER 
S  SCHRAML 
G  SILSBY 
B  SORENSEN 
R  SUMMERS 
W  WALTERS 
AMSRL-WM-TD, 

S  CHOU 
A  M  DIETRICH 
J  M  BOTELER 
D  DANDEKAR 
K  FRANK 
M  RAFTENBERG 
A  RAJENDRAN 
M  SCHEIDLER 
S  SCHOENFELD 
S  SEGLETES  (5  CP) 
T  WEERASOORIYA 
AMSRL-WM-WD, 

J  POWELL 
A  PRAKASH 
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1 


1 


ORGANIZATION 

AERONAUTICAL  &  MARITIME 
RESEARCH  LABORATORY 
ATTN  N  BURMAN 
R  WOODWARD 
S  CIMPOERU 
DPAUL 
PO  BOX  4331 
MELBOURNE  VIC  3001 
AUSTRALIA 

ABTEILUNG  FUER  PHYSHCALISCHE 
CHEMIE 

MONTANUNIVERSITAET 
ATTN  E  KOENIGSBERGER 
A  8700  LEOBEN 
AUSTRIA 

PRB  S  A 

ATTN  M  VANSNICK 

AVENUE  DE  TERVUEREN  168  BTE  7 

BRUSSELS  B  1150 

BELGIUM 

ROYAL  MILITARY  ACADEMY 
ATTN  E  CELENS 
RENAISSANCE  AVE  30 
B1040  BRUSSELS 
BELGIUM 

BULGARIAN  ACADEMY  OF  SCIENCES 
SPACE  RESEARCH  INSTITUTE 
ATTN  VALENTIN  GOSPODINOV 
1000  SOHA  PO  BOX  799 
BULGARIA 

CANADIAN  ARSENALS  LTD 
ATTN  P  PELLETIER 
5  MONTEE  DES  ARSENAUX 
VILLIE  DE  GRADEUR  PQ  J5Z2 
CANADA 

DEFENCE  RSCH  ESTAB  SUFHELD 
ATTN  D  MACKAY 

RALSTON  ALBERTA  TOJ  2NO  RALSTON 
CANADA 


NO.  OF 

COPIES  ORGANIZATION 

1  DEFENCE  RSCH  ESTAB  SUFFIELD 
ATTN  CHRIS  WEICKERT 
BOX  4000  MEDICINE  HAT 
ALBERTA  TIA  8K6 
CANADA 

1  DEFENCE  RSCH  ESTAB  VALCARTIER 
ARMAMENTS  DIVISION 
ATTN  R  DELAGRAVE 
2459  PIE  XI  BLVD  N 
PO  BOX  8800 

CORCELETTE  QUEBEC  GOA  IRO 
CANADA 

1  UNIVERSITY  OF  GUELPH 
PHYSICS  DEPT 
ATTN  C  G  GRAY 
GUELPH  ONTARIO 
NIG  2W1 
CANADA 

1  CEA 

ATTN  ROGER  CHERET 
CEDEX  15 

313  33  RUE  DE  LA  FEDERATION 

PARIS  75752 

FRANCE 

1  CEA 

CISI  BRANCH 
ATTN  PATRICK  DAVID 
CENTRE  DE  SACLAY  BP  28 
GIF  SUR  YVETTE  91192 
FRANCE 

1  CEA/CESTA 

ATTN  ALAIN  GEILLE 
BOX2LE  BARP  33114 
FRANCE 

6  CENTRE  D’ETUDES  DE  GRAMAT 
ATTN  SOLVE  GERARD 
CHRISTIAN  LOUPIAS 
PASCALE  OUTREBON 
J  CAGNOUX 
C  GALLIC 
JTRANCHET 
GRAMAT  46500 
FRANCE 
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2 


3 


6 
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1 
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CENTRE  D’ETUDES  DE  LMEIL-VALENTON 
ATTN  CHRISTIAN  AUSSOURD 
JEAN-CLAUDE  BOZIER 
SAINT  GEORGES  CEDEX 
VILLENEUVE  94195 
FRANCE 

CENTRE  D’ETUDES  DE  VAUJOURS 

ATTN  PLOTARD  JEAN-PAUL 

ERIC  BOTTET 

TAT  SIHN  VONG 

BOITE  POSTALE  NO  7 

COUNTRY  77I8I 

FRANCE 

CENTRE  DE  RECHERCHES 
ET  D’ETUDES  D’ARCUEIL 
ATTN  D  BOUVART 
C  COTTENNOT 
S  JONNEAUX 
H  ORSINI 
S  SERROR 
F  TARDIVAL 

16  BIS  AVENUE  PRIEUR  DE 

LA  COTE  D’OR 

F941 14  ARCUEIL  CEDEX 

FRANCE 

DAT  ETBS  CETAM 

ATTN  CLAUDE  ALTMAYER 

ROUTE  DE  GUERRY  BOURGES 

18015 

FRANCE 

ETBS  DSTI 
ATTN  P  BARNIER 
ROUTE  DE  GUERAY 
BOITE  POSTALE  712 
18015  BOURGES  CEDEX 
FRANCE 

FRENCH  GERMAN  RESEARCH  INST 
ATTN  CHANTERET  P-Y 
CEDEX  12  RUE  DE  I’lNDUSTRIE 
BP  301 

F68301  SAINT-LOUIS 
FRANCE 


5  FRENCH  GERMAN  RESEARCH  INST 
ATTN  HANS-JURGEN  ERNST 
FRANCIS  JAMET 
PASCALE  LEHMANN 
KHOOG 
HLERR 

CEDEX  5  5  RUE  DU  GENERAL 
CASSAGNOU 
SAINT  LOUIS  68301 
FRANCE 

1  LABORATOIRE  DE  TECHNOLOGIE  DES 
SURFACES 

ECOLE  CENTRALE  DE  LYON 
ATTN  VINET  P 
BP  163 

69131  ECULLY  CEDEX 
FRANCE 

1  BATTELLE  INGENIEUTECHNIK  GMBH 
ATTN  W  FUCHE 
DUESSELDORFFER  STR  9 
ESCHBORN  D  65760 
GERMANY 

1  CONDAT 

ATTN  J  KIERMEIR 
MAXIMILIANSTR  28 
8069  SCHEYERN  FERNHAG 
GERMANY 

1  DEUTSCHE  AEROSPACE  AG 
ATTN  MANFRED  HELD 
POSTFACH  13  40 
D  86523  SCHROBENHAUSEN 
GERMANY 

1  DIEHL  GBMH  AND  CO 
ATTN  M  SCHDLDKNECHT 
nSCHBACHSTRASSE  16 
D  90552  ROTBENBACH  AD  PEGNITZ 
GERMANY 
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4  ERNST  MACH  INSTITUT 
ATTN  VOLKER  HOHLER 
E  SCHMOLINSKE 
E  SCHNEIDER 
KTHOMA 
ECKERSTRASSE  4 
D-7800  FREIBURG  I  BR  791  4 
GERMANY 

1  EUROPEAN  SPACE  OPERATIONS  CENTRE 
ATTN  WALTER  FLURY 
ROBERT-BOSCH-STRASSE  5 

64293  DARMSTADT 
GERMANY 

3  FRAUNHOFER  INSTITUT  FUER 
KURZZEITDYNAMIK 
ERNST  MACH  INSTITUT 
ATTN  H  ROTHENHAEUSLER 
HSENF 

E  STRASSBURGER 
HAUPTSTRASSE  18 
D79576  WEIL  AM  RHEIN 
GERMANY 

3  FRENCH  GERMAN  RESEARCH  INST 
ATTN  HARTMUTH  F  LEHR 
ROLF  HUNKLER 
ERICH  WOLLMANN 
POSTFACH  1260 
WEIL  AM  RHEIN  D-79574 
GERMANY 

2  lABG 

ATTN  M  BORRMANN 
H  G  DORSCH 
EINSTEINSTRASSE  20 
D  8012  OTTOBRUN  B  MUENCHEN 
GERMANY 

1  INGENIEURBURO  DEISENROTH 
AUF  DE  HARDT  33  35 
D5204  LOHMAR  1 
GERMANY 


NO.  OF 

COPIES  ORGANIZATION 


3  TU  CHEMNITZ-ZWICKAU 
ATTN  I  FABER 
L  KRUEGER 
LOTHAR  MEYER 

FAKULTAET  FUER  MASCHINENBAU  U. 

VERFAHRENSTECHNIK 
SCHEFFELSTRASSE  110 
09120  CHEMNITZ 
GERMANY 

1  TU  MUNCHEN 
ATTN  E IGENBERGS 
ARCISSTRASSE  21 
8000  MUNCHEN  2 
GERMANY 

2  UNIVERSITAT  PADERBORN 
FACHBEREICH  PHYSIK 
ATTN  O  SCHULTE 

W  B  HOLZAPFEL 
D  33095  PADERBORN 
GERMANY 

1  BHABHA  ATOMIC  RESEARCH  CENTRE 
HIGH  PRESSURE  PHYSICS  DIVISION 
ATTN  N  SURESH 
TROMBAY  BOMBAY  400  085 
INDIA 

1  NATIONAL  OEOPHYSICAL  RESEARCH  INSmUIE 
ATTN  G.  PARTHASARATHY 
HYDERABAD-500  007  (A.  P.) 

INDIA 

5  RAFAEL  BALLISTICS  CENTER 
ATTN  EREZ  DEKEL 
YEHUDA  PARTOM 
G  ROSENBERG 
Z  ROSENBERG 
Y  YESHURUN 
PO  BOX  2250 
HAIFA  31021 
ISRAEL 

1  TECHNION  INST  OF  TECH 
FACULTY  OF  MECH  ENGNG 
ATTN  SOL  BODNER 
TECHNION  CITY 
HAIFA  32000 
ISRAEL 
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I  IHI  RESEARCH  INSTITUTE 
STRUCTURE  &  STRENGTH 
ATTN:  TADASHI  SHIBUE 
1-15,  TOYOSU  3 
KOTO,  TOKYO  135 
JAPAN 

1  ESTEC  CS 

ATTN  DOUGLAS  CASWELL 
BOX  200  NOORDWDK 
2200  AG 
NETHERLANDS 

2  EUROPEAN  SPACE  AGENCY  ESTEC 
ATTN  LUCY  BERTHOUD 
MICHEL  LAMBERT 

POSTBUS  BOX  299  NOORDWDK 

NL2200  AG 

NETHERLANDS 

4  PRINS  MAURTTS  LABORATORY 
ATTN  H  J  REITSMA 
EDWARD  VAN  RIET 
H  PASMAN 
R  YSSELSTEIN 
TNO  BOX  45 
RUSWUK  2280AA 
NETHERLANDS 

1  ROYAL  NETHERLANDS  ARMY 
ATTN  J  HOENEVELD 

V  D  BURCHLAAN  31 
PO  BOX  90822 

2509  LS  THE  HAGUE 
NETHERLANDS 

4  HIGH  ENERGY  DENSITY  RESEARCH  CTR 
ATTN  VLADIMIR  E  FORTOV 
GENADH  I  KANEL 

V  A  SKVORTSOV 
O  YU  VOJOBIEV 
IZHORSKAJA  STR  13/19 
MOSCOW  127412 
RUSSIAN  REPUBLIC 

1  INSTITUTE  OF  CHEMICAL  PHYSICS 
ATTN  A  YU  DOLGOBORODOV 
KOSYGIN  ST  4  V  334 
MOSCOW 

RUSSIAN  REPUBLIC 


3  INSTITUTE  OF  CHEMICAL  PHYSICS 
RUSSIAN  ACADEMY  OF  SCIENCES 
ATTN  A  M  MOLODETS 
S  V  RAZORENOV 
A  V  UTKIN 

142432  CHERNOGOLOVKA 
MOSCOW  REGION 
RUSSIAN  REPUBLIC 

3  rNSHTUTE  OF  MECH  ENGINEERING  PROBLHVIS 

ATTN  V  BULATOV 
D  INDEITSEV 
Y  MESCHERYAKOV 
BOLSHOY,  61,  V.O. 

ST  PETERSBURG  199178 
RUSSIAN  REPUBLIC 

1  INSTITUTE  OF  MINEROLOG  Y  &  PETROGRAPHY 
ATTN  V  A  DREBUSHCHAK 
UNIVERSITETSKI  PROSPEKT,  3 

630090  NOVOSIBIRSK 
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